/*
 *  MWPrint.cpp
 *  T3nsors
 *
 *  Created by Michael Barriault on 10-06-10.
 *  Copyright 2010 MikBarr Studios. All rights reserved.
 *
 */

#define HAVE_HDF

#include "MW.h"
#include <cmath>
#include <fstream>
#include <sdf.h>
using namespace T3;

const real zero = 1e-10;

void MW::Print(real t, Set M) {
#define Phi M.S(0)
#define Psi M.S(1)
#define E M.V(0)
#define B M.V(1)

//	Print(t, Phi);
//	Print(t, Psi);
	Print(t, E);
	Print(t, B);
	real B0 = 1;
	real E0 = 1;
	Scalar EB = E*B/(B0*B0);
	EB.id = "EB";
	Print(t, EB);
	Scalar U = (E*E+B*B)/(E0*E0+B0*B0);
	U.id = "U";
	Print(t, U);
	Scalar divE = div(E);
	Scalar rho = divE/(4*M_PI*E0);
	rho.id = "rho";
	Print(t,rho);
//	Print(t, div(E));
	Print(t, div(B));
	Vector EtB = E%B;
	Scalar EcE = E*curl(E);
	Scalar BcB = B*curl(B);
	Vector J(O,3,"J");
	PFORO FORa J(a)[o] = divE[o]*EtB(a)[o]+(BcB[o]-EcE[o])*B(a)[o];
    Scalar B2 = B*B;
    real B02 = pow(zero*B0,2);
    PFORO FORa if ( B2[o] >  B02 ) J(a)[o] /= 4*M_PI*B2[o]; else J(a)[o] = 0;
    PFORO FORa J(a)[o] *= B0/(E0*E0);
	Print(t, J);
}

void MW::PrintXYZ(real t, Scalar C) {
	real *C2; C2 = &C[0];
	int O2[3];
	FORa O2[a] = O(a);
	real *xyz = new real[O(0)+O(1)+O(2)];
	FOR(i,O(0)) xyz[i] = x(0)(i,0,0);
	FOR(j,O(1)) xyz[O(0)+j] = x(1)(0,j,0);
	FOR(k,O(2)) xyz[O(0)+O(1)+k] = x(2)(0,0,k);
	std::string pid = id + "/0" + C.id;
	gft_out_full((char*)pid.c_str(), t, O2, (char*)"x|y|z", 3, xyz, C2);
	delete xyz;
	return;
}

void MW::PrintL2(real t, Scalar C) {
    real C2[1];
    C2[0] = l2(C);
    int O2[1];
    O2[0] = 1;
    real T[1];
    T[0] = t;
    std::string pid = id + "/l2" + C.id;
    gft_out_full((char*)pid.c_str(), t, O2, (char*)"t", 1, T, C2);
    return;
}

void MW::PrintXY(real t, Scalar C, int k) {
    real *C2;
    C2 = new real[O(0)*O(1)];
    PFOR(i,O(0)) FOR(j,O(1)) C2[i+O(0)*j] = C(i,j,k);
    int O2[2];
    FOR(a,2) O2[a] = O(a);
    real *xy = new real[O(0)+O(1)];
    FOR(i,O(0)) xy[i] = x(0)(i,0,0);
    FOR(j,O(1)) xy[O(0)+j] = x(1)(0,j,0);
    std::string pid = id + "/0" + C.id + "z" + itos((int)x(2)(0,0,k));
    gft_out_full((char*)pid.c_str(), t, O2, (char*)"x|y", 2, xy, C2);
    delete C2;
    delete xy;
    return;
}

void MW::PrintXZ(real t, Scalar C, int j) {
    real *C2;
    C2 = new real[O(0)*O(2)];
    PFOR(i,O(0)) FOR(k,O(2)) C2[i+O(0)*k] = C(i,j,k);
    int O2[2];
    O2[0] = O(0); O2[1] = O(2);
    real *xz = new real[O(0)+O(2)];
    FOR(i,O(0)) xz[i] = x(0)(i,0,0);
    FOR(k,O(2)) xz[O(0)+k] = x(2)(0,0,k);
    std::string pid = id + "/0" + C.id + "y" + itos((int)x(1)(0,j,0));
    gft_out_full((char*)pid.c_str(), t, O2, (char*)"x|z", 2, xz, C2);
    delete C2;
    delete xz;
    return;
}

void MW::PrintYZ(real t, Scalar C, int i) {
    real *C2;
    C2 = new real[O(1)*O(2)];
    PFOR(j,O(1)) FOR(k,O(2)) C2[j+O(1)*k] = C(i,j,k);
    int O2[2];
    O2[0] = O(1); O2[1] = O(2);
    real *yz = new real[O(1)+O(2)];
    FOR(j,O(1)) yz[j] = x(1)(0,j,0);
    FOR(k,O(2)) yz[O(1)+k] = x(2)(0,0,k);
    std::string pid = id + "/0" + C.id + "x" + itos((int)x(0)(i,0,0));
    gft_out_full((char*)pid.c_str(), t, O2, (char*)"y|z", 2, yz, C2);
    delete C2;
    delete yz;
    return;
}

void MW::Print(real t, Vector V) {
	FOR(n,V.N) Print(t, V(n));
	return;
}

void MW::Print(real t, Scalar C) {
    PrintYZ(t, C, O(0)/2);
//    PrintXZ(t, C, O(1)/2);
//    PrintXY(t, C, O(2)/2);
//    PrintXY(t, C, O(2)-1);
    PrintL2(t, C);
}
